
predicted <- c(
  "Belgium", "Botswana", "Burkina Faso", "Czech Republic", "Dominican Republic",
  "Estonia", "France", "Germany", "Greece", "Hungary", "Indonesia", "Japan",
  "Lesotho", "Liberia", "Jamaica", "Mongolia", "Madagascar", "Malawi", "Mozambique",
  "Mexico", "Namibia", "Nigeria", "Nicaragua", "Paraguay", "Poland", "Russia",
  "Slovakia", "Philippines", "Peru", "Switzerland", "South Korea", "Taiwan",
  "Ukraine", "Uganda", "Zimbabwe", "Uruguay", "Denmark", "Netherlands", "Senegal", "Slovenia", "Bulgaria")

opposite <- c(
  "Austria", "Belize", "Benin", "Bolivia", "Brazil", "Cambodia", "Cape Verde",
  "Costa Rica", "Colombia", "Chile", "Cyprus", "Croatia", "Guatemala", "Ireland",
  "Honduras", "Haiti", "Israel", "Italy", "Luxembourg", "Mali", "Portugal",
  "Spain", "Vietnam", "Venezuela")

non_u <- c(
  "Argentina", "Finland", "El Salvador", "Ecuador", "Guyana", "Hong Kong",
  "Kenya", "Malaysia", "Panama", "Sweden", "South Africa", "Suriname",
  "Turkey", "Thailand", "Trinidad & Tobago", "Zambia", "United Kingdom", "Tanzania", "Norway", "Ghana")

df_mu <- read_dta('Data/Reanalysis/jopdataFINAL.dta')

df_predicted <- df_mu %>% filter(newname %in% predicted)

nb.cols <- 41
mycolors <- colorRampPalette(brewer.pal(8, "Set2"))(nb.cols)

m1_pre <- glm(pty_close2cat ~ 0 + distance_pct * newname + I(distance_pct^2) * newname + female + age + 
            urban + high_ed + c1 + c2 + c3 + c4 + c5 + c6 + c7 + c8 + c9 + c10 + c11 + c12 + c13 + 
            c14 + c15 + c16 + c17 + c18 + c19 + c20 + c21 + c22 + c23 + c24 + c25 + c26 + c27 + 
            c28 + c29 + c30 + c31 + c32 + c33 + c34 + c35 + c36 + c37 + c38 + c39 + c40 + c41 + 
            c42 + c43 + c44 + c45 + c46 + c47 + c48 + c49 + c50 + c51 + c52 + c53 + c54 + c55 + 
            c56 + c57 + c58 + c59 + c60 + c61 + c62 + c63 + c64 + c65 + c66 + c67 + c68 + c69 + 
            c70 + c71 + c72 + c73 + c74 + c75 + c76 + c77 + c78 + c79 + c80 + c81 + c82 + c83 + 
            c84 + c85 + c86 + c87 + c88, 
          data = df_predicted, weights = final_weight, family = binomial(link = "logit"))

predicted_plot <- sjPlot::plot_model(m1_pre, type = 'pred', terms = c('distance_pct', 'newname'), ci.lvl = NA)+
  theme_bw()+
  labs(x = 'Distance to/from election',
       y = 'Predictable Probability of Partisanship',
       title = "",
       color = 'Country')+
  scale_color_manual(values=mycolors)+
  theme(legend.position = "none")


df_mu2 <- df_mu %>% filter(newname %in% opposite)

nb.cols <- 25
mycolors <- colorRampPalette(brewer.pal(8, "Set2"))(nb.cols)

m2_opp <- glm(pty_close2cat ~ 0 + distance_pct * newname + I(distance_pct^2) * newname + female + age + 
            urban + high_ed + c1 + c2 + c3 + c4 + c5 + c6 + c7 + c8 + c9 + c10 + c11 + c12 + c13 + 
            c14 + c15 + c16 + c17 + c18 + c19 + c20 + c21 + c22 + c23 + c24 + c25 + c26 + c27 + 
            c28 + c29 + c30 + c31 + c32 + c33 + c34 + c35 + c36 + c37 + c38 + c39 + c40 + c41 + 
            c42 + c43 + c44 + c45 + c46 + c47 + c48 + c49 + c50 + c51 + c52 + c53 + c54 + c55 + 
            c56 + c57 + c58 + c59 + c60 + c61 + c62 + c63 + c64 + c65 + c66 + c67 + c68 + c69 + 
            c70 + c71 + c72 + c73 + c74 + c75 + c76 + c77 + c78 + c79 + c80 + c81 + c82 + c83 + 
            c84 + c85 + c86 + c87 + c88, 
          data = df_mu2, weights = final_weight, family = binomial(link = "logit"))

opposite_plot <- sjPlot::plot_model(m2_opp, type = 'pred', terms = c('distance_pct', 'newname'), ci.lvl = NA)+
  theme_bw()+
  labs(x = 'Distance to/from election',
       y = 'Predictable Probability of Partisanship',
       title = "",
       color = 'Country')+
  scale_color_manual(values=mycolors)+
  theme(legend.position = "none")



df_mu3 <- df_mu %>% filter(newname %in% non_u)

nb.cols <- 21
mycolors <- colorRampPalette(brewer.pal(8, "Set2"))(nb.cols)

m3_non <- glm(pty_close2cat ~ 0 + distance_pct * newname + I(distance_pct^2) * newname + female + age + 
            urban + high_ed + c1 + c2 + c3 + c4 + c5 + c6 + c7 + c8 + c9 + c10 + c11 + c12 + c13 + 
            c14 + c15 + c16 + c17 + c18 + c19 + c20 + c21 + c22 + c23 + c24 + c25 + c26 + c27 + 
            c28 + c29 + c30 + c31 + c32 + c33 + c34 + c35 + c36 + c37 + c38 + c39 + c40 + c41 + 
            c42 + c43 + c44 + c45 + c46 + c47 + c48 + c49 + c50 + c51 + c52 + c53 + c54 + c55 + 
            c56 + c57 + c58 + c59 + c60 + c61 + c62 + c63 + c64 + c65 + c66 + c67 + c68 + c69 + 
            c70 + c71 + c72 + c73 + c74 + c75 + c76 + c77 + c78 + c79 + c80 + c81 + c82 + c83 + 
            c84 + c85 + c86 + c87 + c88, 
          data = df_mu3, weights = final_weight, family = binomial(link = "logit"))

nonu_plot <- sjPlot::plot_model(m3_non, type = 'pred', terms = c('distance_pct', 'newname'), ci.lvl = NA)+
  theme_bw()+
  labs(x = 'Distance to/from election',
       y = 'Predictable Probability of Partisanship',
       title = "",
       color = 'Country')+
  scale_color_manual(values=mycolors)+
  theme(legend.position = "none")


MU_plot <- suppressMessages(ggpubr::ggarrange(predicted_plot, opposite_plot, nonu_plot, 
                             labels = c("A", "B", "C"),
                             ncol = 2, nrow = 2))

suppressMessages(print(MU_plot))

#ggsave(MU_plot, file = 'MU_plot.png', units = 'in', width = 8, height = 8, dpi = 350)
